myhist3D <- function(z,
                     x,
                     y,
                     zlab = NULL,
                     xlab = NULL,
                     ylab = NULL,
                     zlim = NULL,
                     xlim = NULL,
                     ylim = NULL,
                     nlines = 5,
                     cols = NULL,
                     border = NA,
                     main = NULL,
                     cex.main = 1,
                     ntick = NULL,
                     expand = 0.8,
                     theta = 45,
                     phi = 0,
                     r = sqrt(3),
                     d = 1.2,
                     add = FALSE) {
  if (!("plot3D" %in% installed.packages())) {
    warning("Package 'plot3D' required. Installing now.")
    install.packages("plot3D")
    library("plot3D")
  }
  if (!(is.matrix(z))) {
    stop("z: Please supply a matrix for z.")
  }
  if (!(is.vector(x) &
        is.numeric(x) & length(x) == dim(z)[1] + 1)) {
    stop("x: Please supply a numeric vector of the length of the first
         dimension of z plus one.")
  }
  if (!(is.vector(y) &
        is.numeric(y) & length(y) == dim(z)[2] + 1)) {
    stop("y: Please supply a numeric vector of the length of the second
         dimension of z plus one.")
  }
  
  ## Adjust x and y
  x <- x[-length(x)] + abs(x[1] - x[2]) / 2
  y <- y[-length(y)] + abs(y[1] - y[2]) / 2
  
  length.x <- length(x)
  length.y <- length(y)
  nlines <- nlines
  points.x <- x[round(seq(1, length(x), length.out = nlines))]
  points.y <- y[round(seq(1, length(y), length.out = nlines))]
  if (is.null(ntick)) {
    ntick <- nlines
  }
  if (is.null(zlim)) {
    zlim <- range(z)
  }
  if (is.null(xlim)) {
    xlim <- range(x)
  }
  if (is.null(ylim)) {
    ylim <- range(y)
  }
  if (any(is.na(cols))) {
    col <- ramp.col (col = c("gray60", "gray10"),
                     n = 100,
                     alpha = .9)
  } else {
    col <- ramp.col (col = cols, n = 100, alpha = .9)
  }
  
  if (isFALSE(add)) {
    perspbox(
      z = NULL,
      zlim = zlim,
      zlab = zlab,
      x = NULL,
      xlim = xlim,
      xlab = xlab,
      y = NULL,
      ylim = ylim,
      ylab = ylab,
      main = main,
      cex.main = cex.main,
      bty = 'u',
      col.grid = "white",
      col.panel = "gray95",
      col.axis = "gray80",
      lwd.grid = 2,
      lwd.panel = 2,
      ticktype = "detailed",
      ntick = ntick,
      expand = expand,
      theta = theta,
      phi = phi,
      r = r,
      d = d
    )
  }
  
  hist3D (
    x = x,
    y = y,
    z = z,
    col = col,
    border = border,
    add = TRUE,
    plot = TRUE
  )
}